home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / MRQCOF.DEM < prev    next >
Text File  |  1991-04-29  |  3KB  |  110 lines
  1. PROGRAM d14r9(input,output);
  2. (* driver for routine MRQCOF *)
  3. CONST
  4.    npt=100;
  5.    mma=6;
  6.    spread=0.1;
  7. TYPE
  8.    glndata = ARRAY [1..npt] OF real;
  9.    glmma = ARRAY [1..mma] OF real;
  10.    glnparam=glmma;
  11.    gllista = ARRAY [1..mma] OF integer;
  12.    glnalbynal = ARRAY [1..mma,1..mma] OF real;
  13. VAR
  14.    glinext,glinextp : integer;
  15.    glma : ARRAY [1..55] OF real;
  16.    gliset : integer;
  17.    glgset : real;
  18.    chisq : real;
  19.    i,j,idum,mfit : integer;
  20.    x,y,sig : glndata;
  21.    a,beta,gues : glmma;
  22.    lista : gllista;
  23.    covar,alpha : glnalbynal;
  24.  
  25. (*$I MODFILE.PAS *)
  26. (*$I RAN3.PAS *)
  27.  
  28. (*$I GASDEV.PAS *)
  29.  
  30. PROCEDURE funcs(x: real; a: glnparam; VAR y: real;
  31.        VAR dyda: glnparam; na: integer);
  32. (* Programs using routine FGAUSS must define the type
  33. TYPE
  34.    glnparam = ARRAY [1..na] OF real;
  35. where na is the number of parameters. *)
  36. VAR
  37.    i,ii : integer;
  38.    fac,ex,arg : real;
  39. BEGIN
  40.    y := 0.0;
  41.    FOR ii := 1 to (na DIV 3) DO BEGIN
  42.       i := 3*ii-2;
  43.       arg := (x-a[i+1])/a[i+2];
  44.       ex := exp(-sqr(arg));
  45.       fac := a[i]*ex*2.0*arg;
  46.       y := y+a[i]*ex;
  47.       dyda[i] := ex;
  48.       dyda[i+1] := fac/a[i+2];
  49.       dyda[i+2] := fac*arg/a[i+2]
  50.    END
  51. END;
  52.  
  53. (*$I MRQCOF.PAS *)
  54.  
  55. BEGIN
  56.    gliset := 0;
  57.    a[1] := 5.0; a[2] := 2.0; a[3] := 3.0;
  58.    a[4] := 2.0; a[5] := 5.0; a[6] := 3.0;
  59.    gues[1] := 4.9; gues[2] := 2.1; gues[3] := 2.9;
  60.    gues[4] := 2.1; gues[5] := 4.9; gues[6] := 3.1;
  61.    idum := -911;
  62. (* first try sum of two gaussians *)
  63.    FOR i := 1 to 100 DO BEGIN
  64.       x[i] := 0.1*i;
  65.       y[i] := 0.0;
  66.       y[i] := y[i]+a[1]*exp(-sqr((x[i]-a[2])/a[3]));
  67.       y[i] := y[i]+a[4]*exp(-sqr((x[i]-a[5])/a[6]));
  68.       y[i] := y[i]*(1.0+spread*gasdev(idum));
  69.       sig[i] := spread*y[i]
  70.    END;
  71.    mfit := mma;
  72.    FOR i := 1 to mfit DO BEGIN
  73.       lista[i] := i
  74.    END;
  75.    FOR i := 1 to mma DO BEGIN
  76.       a[i] := gues[i]
  77.    END;
  78.    mrqcof(x,y,sig,npt,a,mma,lista,mfit,alpha,beta,mma,chisq);
  79.    writeln;
  80.    writeln('matrix alpha');
  81.    FOR i := 1 to mma DO BEGIN
  82.       FOR j := 1 to mma DO write(alpha[i,j]:12:4);
  83.       writeln
  84.    END;
  85.    writeln('vector beta');
  86.    FOR i := 1 to mma DO write(beta[i]:12:4);
  87.    writeln;
  88.    writeln('chi-squared:',chisq:12:4);
  89.    writeln;
  90. (* next fix one line and improve the other *)
  91.    FOR i := 1 to 3 DO BEGIN
  92.       lista[i] := i+3
  93.    END;
  94.    mfit := 3;
  95.    FOR i := 1 to mma DO BEGIN
  96.       a[i] := gues[i]
  97.    END;
  98.    mrqcof(x,y,sig,npt,a,mma,lista,mfit,alpha,beta,mma,chisq);
  99.    writeln('matrix alpha');
  100.    FOR i := 1 to mfit DO BEGIN
  101.       FOR j := 1 to mfit DO write(alpha[i,j]:12:4);
  102.       writeln
  103.    END;
  104.    writeln('vector beta');
  105.    FOR i := 1 to mfit DO write(beta[i]:12:4);
  106.    writeln;
  107.    writeln('chi-squared:',chisq:12:4);
  108.    writeln
  109. END.
  110.